
***************************************************
* Figure 3
***************************************************


use "$workdata/parents", clear /* see $data_do/parents.do */

g b_month=month(foed_dag)
g b_year=year(foed_dag)

*****************************************************************

collapse wage (sum) e, by(pnr age_q month year b_month b_year parent) 


capture drop age_month_q_y
g age_month_q_y=floor(age_q)
tostring age_month_q_y, replace
capture drop age_month_q_m
g age_month_q_m=(round((age_q-floor(age_q))*12)+1)/100
tostring age_month_q_m, replace
capture drop age_month_q
egen age_month_q=concat(age_month_q_y age_month_q_m)
destring age_month_q, replace
local resp="age_month_q>=24.12 & age_month_q<=25.04"
g resp=`resp'
g d=age_q>=25
g age_q0=age_q-25
g age_q0_2= age_q0^2
g age_q0_3= age_q0^3
g dXage_q0=d*age_q0
g dXage_q0_2=d*age_q0_2
g dXage_q0_3=d*age_q0_3

g np=parent==0
g treat=d

g age_q0Xnp=age_q0*np
g age_q0_2Xnp=age_q0_2*np
g age_q0_3Xnp=age_q0_3*np
g dXage_q0Xnp=d*age_q0*np
g dXage_q0_2Xnp=d*age_q0_2*np
g dXage_q0_3Xnp=d*age_q0_3*np
g age_q0XtreatXnp=age_q0*treat*np
g age_q0_2XtreatXnp=age_q0_2*treat*np
g age_q0_3XtreatXnp=age_q0_3*treat*np
g treatXnp=treat==1 & np==1

*********************************************************************

eststo clear
local model=0
foreach y of varlist wage {
forvalues s=2/2  {

preserve
keep if age_q>=21 & age_q<29


reg `y' i.year#i.month i.b_year#i.b_month [aw=e] if parent==0
predict res_parent0, r
replace `y'=res_parent0 if parent==0


reg `y' i.year#i.month i.b_year#i.b_month [aw=e] if parent==1
predict res_parent1, r
replace `y'=res_parent1 if parent==1



local spec_num=`s'
if `y'==wage local title="Employment (residuals)"

local min_wage=-0.06
local max_wage=0.02
local min_wage_p=-0.02
local max_wage_p=0.06


if `spec_num'==2 local spec="d age_q0 age_q0_2          dXage_q0 dXage_q0_2"
if `spec_num'==2 local spec_name="2nd order polynomial spline"
local model =`model'+1
reg `y' `spec' [aw=e] if resp==0 & parent==0, vce(cluster pnr)
predict yhat 
local beta_d=round(_b[d], .001)*100

reg `y' `spec' [aw=e] if  parent==1, vce(cluster pnr)
predict yhat_p 
local beta_d_p=round(_b[d], .001)*100

local spec_num=`s'
if `y'==wage local title="Employment"
if `y'==wage local pos=1  /*11*/
if `y'==wage local place="neast"  /*"nwest"*/

collapse `y' yhat yhat_p (sem) se=`y'  [aw=e] if age_q>=21 & age_q<29, by(age_q parent)
g `y'_min=`y'-se*1.96 if parent==0
g `y'_max=`y'+se*1.96 if parent==0
g `y'_min_p=`y'-se*1.96 if parent==1
g `y'_max_p=`y'+se*1.96 if parent==1

capture drop age_month_q_y
g age_month_q_y=floor(age_q)
tostring age_month_q_y, replace
capture drop age_month_q_m
g age_month_q_m=(round((age_q-floor(age_q))*12)+1)/100
tostring age_month_q_m, replace
capture drop age_month_q
egen age_month_q=concat(age_month_q_y age_month_q_m)
destring age_month_q, replace
local resp="age_month_q>=24.12 & age_month_q<=25.04"
g resp=`resp'


twoway (rarea `y'_min `y'_max age_q if age_month_q<24.12 & parent==0, color(gs13%80) yaxis(1) ) ///
	   (rarea `y'_min `y'_max age_q if age_month_q>25.04 & parent==0, color(gs13%80) yaxis(1)) ///
	   (scatter `y' age_q if age_q< 25 & parent==0                        /*& `y'>`min_`y'' & `y'<`max_`y''*/, msymbol(o) mlcolor(gs3) mfcolor(gs16)  yaxis(1))  ///
	   (scatter `y' age_q if age_q>=25 & parent==0                        /*& `y'>`min_`y'' & `y'<`max_`y''*/, msymbol(o) mlcolor(gs3) mfcolor(gs16)  yaxis(1)) ///
	   (scatter `y' age_q if `resp'    & parent==0                                                           , msymbol(o) mlcolor(gs3) mfcolor(gs16)  yaxis(1))     ///
	   (line yhat age_q if age_q< 25   & parent==0                        /*& `y'>`min_`y'' & `y'<`max_`y''*/, lcolor(gs3%90) lpattern(solid) yaxis(1)) ///
   	   (line yhat age_q if age_q>= 25  & parent==0                        /*& `y'>`min_`y'' & `y'<`max_`y''*/, lcolor(gs3%90) lpattern(solid) yaxis(1)) ///
	   (rarea `y'_min_p `y'_max_p age_q if age_month_q<=24.12 & parent==1  & `y'>`min_`y'_p' & `y'<`max_`y'_p' , color(gs13%80) yaxis(2)) ///
	   (rarea `y'_min_p `y'_max_p age_q if age_month_q>=25.01 & parent==1  & `y'>`min_`y'_p' & `y'<`max_`y'_p' , color(gs13%80) yaxis(2) ) ///
	   (scatter `y' age_q if age_q< 25 & parent==1  					  & `y'>`min_`y'_p' & `y'<`max_`y'_p', msymbol(O) mlcolor(gs3) mfcolor(gs16) /*mcolor(gs3%90)*/ yaxis(2) )  ///
	   (scatter `y' age_q if age_q>=25 & parent==1  					  & `y'>`min_`y'_p' & `y'<`max_`y'_p', msymbol(O) mlcolor(gs3) mfcolor(gs16) /*mcolor(gs3%90)*/ yaxis(2)) ///
	   (scatter `y' age_q if `resp'    & parent==1 						  & `y'>`min_`y'_p' & `y'<`max_`y'_p', msymbol(O) mlcolor(gs3) mfcolor(gs16) /*mcolor(gs8%90)*/ yaxis(2)) /// 
	   (line yhat_p age_q if age_q< 25 & parent==1                        & `y'>`min_`y'_p' & `y'<`max_`y'_p', lcolor(gs3) lpattern(dash) yaxis(2)) ///
   	   (line yhat_p age_q if age_q>=25 & parent==1  					  & `y'>`min_`y'_p' & `y'<`max_`y'_p', lcolor(gs3) lpattern(dash) yaxis(2)) ///
	   , subtitle(`title' (residuals), margin(b=4)) ytitle("Employment, Non-Parents", height(5) axis(1)) ylabel(`min_`y''  (0.02)`max_`y''  , axis(1))  yscale(r(-0.07 0.03) titlegap(2) axis(1))   ///
	     ytitle("Employment, Parents", height(5) axis(2)) ylabel(`min_`y'_p'(0.02)`max_`y'_p', axis(2))  yscale(r(-0.03 0.07) titlegap(2) axis(2)) ///
	     xlabel(21(1)29)  xtitle("Age in months", height(5))   xline(24.958333) ///
		 scheme(s1mono)  name(fig3_1, replace) note("Non-Parents. Effect size: `beta_d' pp (Specification: `spec_name')" "Parents. Effect size: `beta_d_p' pp (Specification: `spec_name')" ) ///
		legend(order( 3 "Non-Parents (left axis)" 10 "Parents (right axis)") lcolor(gs6) lwidth(thin) col(1) position(`pos') ring(0) bplace(`place')) 

graph save fig3_1  $figures/fig3_1, replace
graph export $figures/fig3_1.pdf, name(fig3_1) replace
graph export $figures/fig3_1.png, name(fig3_1) replace

sum yhat if age_q<25 & parent==0
local parent0_m=r(mean)
di `parent0_m'
replace yhat=yhat-`parent0_m'
sum yhat_p if age_q<25 & parent==1
local parent1_m=r(mean)
di `parent1_m'
replace yhat_p=yhat_p-`parent1_m'


twoway (line yhat age_q if age_q< 25   & parent==0                        , lcolor(gs3%90) lpattern(solid)) ///
   	   (line yhat age_q if age_q>= 25  & parent==0                        , lcolor(gs3%90) lpattern(solid)) ///
	   (line yhat_p age_q if age_q< 25 & parent==1                        & `y'>`min_`y'_p' & `y'<`max_`y'_p', lcolor(gs3) lpattern(dash) ) ///
   	   (line yhat_p age_q if age_q>=25 & parent==1  					  & `y'>`min_`y'_p' & `y'<`max_`y'_p', lcolor(gs3) lpattern(dash) ) ///
	   , subtitle(`title' (regression lines), margin(b=4)) ytitle("Employment (residuals, centered)", height(5) ) ylabel(-0.04 (0.02) 0.04 )    ///
	     xlabel(21(1)29)  xtitle("Age in months", height(5))   xline(24.958333) ///
		 scheme(s1mono)  name(fig3_2, replace) note("Non-Parents. Effect size: `beta_d' pp (Specification: `spec_name')" "Parents. Effect size: `beta_d_p' pp (Specification: `spec_name')" ) ///
		legend(order( 1 "Non-Parents" 3 "Parents") lcolor(gs6) lwidth(thin) col(1) position(`pos') ring(0) bplace(`place')) 

graph save fig3_2  $figures/fig3_2, replace
graph export $figures/fig3_2.pdf, name(fig3_2) replace
graph export $figures/fig3_2.png, name(fig3_2) replace

graph use $figures/fig3_1
graph use $figures/fig3_2
graph combine ///
      fig3_1 fig3_2 ///
      , cols(2) rows(2) ysize(4) xsize(10) altshrink  scheme(s1mono) ///
	   name(fig3, replace)
graph save fig3  $figures/fig3, replace
graph export $figures/fig3.pdf, name(fig3) replace
graph export $figures/fig3.png, name(fig3) replace

restore
}
}

***************************************************
* end: Figure 3
***************************************************
